//----------------------------------------------------------------------------
// M A R I T I M E  R E S E A R C H  I N S T I T U T E  N E T H E R L A N D S
//----------------------------------------------------------------------------
// Copyright (C) 2008 - MARIN - All rights reserved - http://www.marin.nl
//----------------------------------------------------------------------------
// Program    : mermaid
// Module     : vbm
// File       : Solver.cpp
//----------------------------------------------------------------------------

#ifdef _WINDOWS
#define NOMINMAX
#endif

#include "Solver.h"
#include "FloatUtils.h"
#include "Timer.h"

#include <algorithm>
#include <fstream> // voor printf

#ifdef EWLAPACK
#include <lapackpp.h> // Lapack++ library // see NopDefSolver.cpp
#include <sybfd.h> // for LaSymmBandMatFactorizeIP()
#endif

#define COLUMNS  (unsigned)(Nx1 + 3)
#define ROWS     (unsigned)(Nx2 + 3)

/// constructor
Solver::Solver(int Nx1, int Nx2, REAL tolerance)
: m_Nx1(Nx1)
, m_Nx2(Nx2)
, m_tolerance(tolerance)
// some member variables for the deflation method
{
}

/// destructor
Solver::~Solver()
{
}

#ifdef EWLAPACK
// A Lapack++ routine for solving a linear system with spd matrix.
void Solver::LaLinearSolveIPChol(LaSymmBandMatDouble &A,
                                 LaGenMatDouble &B)
{
    /* This function is copied from the LaLinearSolveIP() function of Lapack++.
     * The original function calculates x from Ax=b with A=LL^T the Cholesky decomposition.
     * As input, the matrix A and right hand side b are needed. In this function, the same
     * system Ax=b is solved, with A=LL^T the Cholesky decomposition. However, the input is
     * now given by the lower triangular matrix L and the righ hand side b. */

    assert(A.size(1) == B.size(0));

    //LaSymmBandMatFactorizeIP(A); // This statement is removed from LaLinearSolveIP().

    char uplo = 'L';
    integer N = A.size(1);
    integer KD = A.subdiags();
    integer M = B.size(1);
    integer lda = A.gdim(0);
    integer ldb = B.gdim(0);
    integer info = 0;

    F77NAME(dpbtrs)(&uplo, &N, &KD, &M, &A(0, 0), &lda, &B(0, 0), &ldb, &info);

    assert(info == 0);
} // end of LaLinearSolveIPChol()
#endif
